Numerical evolution of Brill waves 



David Garfinkletl 

Department of Physics, Oakland University, Rochester, Michigan 48309 

G. Comer Duncani 

Department of Physics and Astronomy, Bowling Green State University, Bowling Green, Ohio 43403 

(February 7, 2008) 



O 
O 
O 

(N 
•*-> 
O 

O 
o 

(N 

(N 
> 

m 
o 

VO 
O 
O 

o 

I ' 



x 



We report a numerical evolution of axisymmetric Brill waves. The numerical algorithm has new 
features, including (i) a method for keeping the metric regular on the axis and (ii) the use of 
coordinates that bring spatial infinity to the edge of the computational grid. The dependence of the 
evolved metric on both the amplitude and shape of the initial data is found. 
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I. INTRODUCTION 

In the last decade, there has been much use of numer- 
ical simulations to study gravitational collapse. These 
studies include critical collapse 0,0, black hole colli- 
sions HHH] and the approach to the singularity (|^,||,^| ■ 
In simulations, spherical or planar symmetry allows the 
finest resolution, while complete lack of symmetry (a full 
3+1 simulation) would allow one to treat a completely 
general situation. An intermediate case is axisymmetry, 
which allows more resolution than a full 3+1 simulation; 
but also permits the study of situations (e.g. prolate 
collapse or black hole collisions) that do not occur in 
spherical symmetry. 

One issue that can be addressed in axisymmetry is 
weak cosmic censorship: whether the singularities formed 
in gravitational collapse are hidden inside black hole 
event horizons. A simulation done by Shapiro and 
Teukolsky Jl(| of the collapse of collisionless matter indi- 
cates that weak cosmic censorship may be violated in the 
collapse of highly prolate objects. The result of reference 
JT(| is that for highly prolate initial configurations, a sin- 
gularity forms before an apparent horizon. This result 
might still be consistent with weak cosmic censorship if 
it is an artifact of the slicing |ll[] or the type of matter 
used. To address the second possibility, one would like 
to know whether the same type of behavior occurs in 
vacuum collapse. 

In ]l^| Abrahams et al examine families of initial data 
for Brill waves, fig] ] These are vacuum, axisymmetric ini- 
tial data at a moment of time symmetry. The authors 
of Jl^] show that by considering sufficiently prolate con- 
figurations, one can find initial data with no apparent 
horizon but with large values of the Riemann invariant 
R abcd R a bcd- They conjecture that such initial data, when 
evolved, will form singularities without apparent hori- 
zons. 

In order to test this conjecture, one would like to evolve 
initial data for highly prolate Brill waves to find the be- 
havior of the evolved spacetime. More generally, one 



would like to know how the collapse process depends on 
both the shape and the strength of the initial data. In 
this paper, we report numerical simulations of the col- 
lapse of Brill waves. We find the dependence of the col- 
lapse on the amplitude and the initial shape of the wave. 
The numerical method is presented in section 2 and the 
results in section 3. Conclusions are given in section 4. 



II. NUMERICAL METHOD 

One difficulty present in numerical simulations of ax- 
isymmetric spacetimes is the existence of the axis where 
the Killing field vanishes. If one chooses coordinates 
adapted to the symmetry, then there is a coordinate sin- 
gularity on the axis. This coordinate singularity has the 
possibility of causing numerical problems: that is, one 
must be very careful in the choice of numerical method to 
ensure that the metric remains smooth on the axis. Most 
axisymmetric simulations use spherical coordinates. In 
addition to a coordinate singularity on the axis, spherical 
coordinates lead to a more severe coordinate singularity 
at the origin. This singularity is either avoided by using 
initial data with a black hole and no origin, Q or treated 
by using elaborate numerical methods to keep the metric 
regular at the origin. |l4| 

We bypass this origin difficulty by the use of cylindrical 
coordinates (z,r,cf>). The spatial metric j a b takes the 
form 



ds 2 = yj 4 [e 2rS (dz 2 + dr 2 ) + r 2 



(1) 



Here, ip and S are functions of z, r and the time t. The 
axis is at r = and is the only coordinate singularity. 

Numerically, the axis is an edge of the computational 
grid, and therefore values for the variables must be given 
on the axis points. However, analytically the axis con- 
sists of interior points, and therefore the only permiss- 
ablc condition to impose is smoothness. For a scalar /, 
smoothness requires that / be an even function of r and 
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therefore that d r f vanish on the axis. This condition 
can be used numerically to set the value of / at the axis 
points. Similarly, for components of tensors in cylindrical 
coordinates, smoothness requires either that the compo- 
nent be an even function of r or that it be an odd function 
of r. Odd functions of r vanish on the axis, while even 
functions have vanishing r derivative on the axis. In ci- 
ther case this information can be used to set the value of 
the variable at the axis points. 

Unfortunately, a difficulty arises because smoothness 
requires that certain quantities be even and also vanish 
on the axis. One such quantity is K r r — K^^. This quan- 
tity must both vanish and have vanishing derivative on 
the axis. However, since we can impose only one bound- 
ary condition at the edge of the computational grid, we 
do not have a way to maintain both conditions. Our so- 
lution to this difficulty is to introduce a new variable: de- 
fine W = (K r r —K^ 0) jr. Then W is odd and smoothness 
requires only that W vanish on the axis. This method is 
the reason that the term in the exponential in equation 
(|l|) is written as 2rS. Smoothness of the spatial metric 
requires that 7^/ '{r 2 j rr ) = 1 + o(r 2 ) and therefore that 
the quantity rS both vanish and have vanishing deriva- 
tive on the axis. This condition is met if S vanishes on 
the axis. In the end, all quantities we deal with either 
vanish on the axis or have vanishing derivative there, but 
not both. If we encounter a quantity / that is order r 2 at 
the axis, we simply define the odd quantity g = f/r and 
rewrite all equations containing / in terms of g. 

There remains the question of how to numerically im- 
plement the appropriate conditions on the axis points. 
For an odd quantity X , the most natural method would 
be to put the first grid point at r = and impose the 
condition X(l) — 0. However, for an even quantity Y, 
the simple condition Y(l) = Y(2) would then impose 
d r Y = not on the axis but at r = A r /2 where A r 
is the grid spacing in r. Instead we use the method of 
"ghost zones": the first grid point is at r — —A r /2. The 
axis is then halfway between gridpoints 1 and 2. For an 
odd quantity X, we impose the condition X(l) = —X(2) 
while for an even quantity Y, we use Y(l) = Y(2). In 
each case, the appropriate condition is satisfied on the 
axis. 

We note that our method is not the only way to keep 
the axis stable. The "cartoon" method of reference ]lq ] 
begins with a cartesian 3+1 code and operates it in a 
thin slab with boundary conditions at the faces of the 
slab given by the axisymmctry of the solution. This car- 
toon method is effective (see reference |l6) for another 
implementation). However, since a 3+1 code uses more 
variables than an axisymmetric code, for a given axisym- 
metric problem the cartoon method uses more computer 
memory than our method. 

We now turn to the method of evolution. The spatial 
metric jab is evolved using the ADM equation 

dtlab = -2aK ab + C(3-f ab (2) 

where K ab is the extrinsic curvature, a is the lapse and (3 a 



is the shift. From the form of the metric in equation (|l|) it 
is clear that we have imposed the conditions j rz — and 
j rr = j zz . In order that these conditions be preserved 
by the evolution in equation (O), the shift must satisfy 



d r f3 z + d z p r = 2aK z 



(3) 



d z (3 z - d r p r = aU (4) 
where U = K z z — K r r . This gives rise to the equations 
d r d r p r + d z d z (3 r = 2d z (aK z r ) - d r (aU) (5) 

d r d r (3 z + d z d z (3 z = 2d r (aK z r ) + d z {aU) (6) 

Equation (g) yields the following evolution equation 
for S 

d t S = -aW + /3 z d z S + p r d r S + [i r S/r + d r {(3 r /r) (7) 

where W = (K r r — K^^)/r. Rather than evolve ip, we 
solve for it using the Hamiltonian constraint. 

We choose maximal slicing {K = 0) and evolve the 
extrinsic curvature using the ADM equation 



d t K a b = -D a D b a + aR a b + CpK a t 



(8) 



where R a b is the Ricci tensor of the spatial metric. Since 
K = 0, the only independent components of the extrinsic 
curvature are K z r , U and W . These evolve as follows: 

d t K z r = ^e~ 2rS [(S + 2i\)- x d r ^ + rd r S) d z a 
—d z d r a + (2ip~ 1 d z ip + rd z S) d r a + aR zr ] 

+(3 z 8 z K z r + ffd r K z r + hj (d r p z - d z (3 r ) (9) 

d t U = t/j- A e- 2rS [(4t/;- 1 ^V + 2rd z S) d z a 
- (2S + A^dri) + 2rd r S) d r a 
+d r d r a — d z d z a — aR a ] 
+(3 z d z U + (3 r d r U + 2K z r (d z (3 r - d r (3 z ) (10) 

d t W = tp- 4 e- 2rS [-d r {r- l d r a) - d z Sd z a 
+ (s/r + d r S + 4(r^r 1 d r v) d r a] 
+aR b + (3 z d z W + (3 r d r W + W/3 r /r 

+ {K Z r/r) (8 r p z - d z (3 r ) (11) 

Here we have R a = R rr — R zz and R b = (R r r — R^^/r. 

The evolution must preserve the condition K = 0, 
which implies that the lapse must satisfy D a D a a = 
aK a b K b a . This equation becomes 



aip e 



6 2rS 



r 1 d r (ri/j 2 d r a) + d z (■0 2 d z a) 
\{U 2 + r 2 W 2 + UrW) + 2{K z r f 



(12) 
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The Hamiltonian constraint, R—K a t,K a — becomes 
the following equation for the conformal factor ip. 

(d r d r + r~~ 1 d r + d z d z )ip 



(d r d r + d z d z ){rS) 



+^e 2rS ( ^(U* + r z W' z + UrW) + (K z r )' 



(13) 



Our set of variables is then (S, ip, K z r , U, W, a, pf ', f3 z ). 
Of these variables, S, K z r , W and [3 r are odd functions of 
r, while the rest are even functions. These variables are 
evolved as follows: S and the extrinsic curvature vari- 
ables are evolved using equations ( ^ |Io| , |Tl"| ) . At each 
time step, the elliptic equations (f5|j6|,12|,|l3|) are solved for 
the shift, lapse and conformal factor. 

To begin the evolution, we need initial data satisfying 
the constraint equations. Initial data for Brill waves is 
a moment of time symmetry, so K a t — and therefore 
the momentum constraint is automatically satisfied. The 
variable S can be freely specified (subject to smoothness 
on the axis and asymptotic flatness at infinity). Our 
choice for S is 



ar exp 



(14) 



where a, <r r and a z are constants. Here, a is the ampli- 
tude of the wave and a r and a z are widths in the r and z 
directions respectively. Given S, equation ( |l3| ) is solved 
for ip. 

There remains the question of the boundary conditions 
to apply at the outer edge of the computational grid. A 
natural method would be to put the outer edge of the 
computational grid at some large distance and to impose 
some outgoing wave condition on the evolution equations 
and a Robin boundary condition on the elliptic equa- 
tions. However, the issue of appropriate boundary con- 
ditions for a mixed hyperbolic-elliptic set of equations is 
quite complicated. This issue becomes even more dificult 
in cylindrical coordinates than in spherical coordinates 
since the asymptotic behavior of the variables looks more 
complicated in cylindrical coordinates. Our attempts to 
impose a boundary condition of this sort led to numer- 
ical instability. Instead, we decided to use a different 
approach. We begin by noting that a coordinate transfor- 
mation can bring spatial infinity to the edge of the com- 
putational grid. We introduce new coordinates (z, r) de- 
fined by z = tan z and r = tan f. We then place the edges 
of the computational grid at z = it/2 and at r = tt/2. 
These regions correspond to spatial infinity. Though we 
use new coordinates, we retain the old metric and ex- 
trinsic curvature variables, with the exception that we 
introduce the quantities S = S/ cos r and W = W/ cos f. 
Thus our set of variables is (S, ip, K z r , U, W, a, ff , (3 Z ) 
and the set of equations used to evolve these variables 
is equations (@,|Jl^Jl|,|, HH© with the substitutions 



r — > tanf, d r — > cos fdf, d z — > cos zd z , S — > cosrS and 
W — * cosfW. 

Since the edge of the grid is at spatial infinity, the 
choice of boundary condition is dictated by asymptotic 
flatness: ip and a must be 1 at the outer boundary, and 
all the other variables must vanish there. Though the co- 
ordinate transformation is singular, this does not lead to 
a singularity in the evolution equations. In fact, just the 
opposite is true: in the new coordinates, the right hand 
sides of the evolution equations approach zero at spatial 
infinity (as they must, since the variables are unchanging 
there). 

The advantages of this spatial infinity boundary condi- 
tion are stability and consistency with the field equations. 
However, there are disadvantages as well. The change of 
variables changes the form of the elliptic operators that 
need to be inverted to solve the elliptic equations. This 
slows down the process of solving the elliptic equations 
and therefore slows down the code. A more fundamental 
difficulty has to do with waves produced in the collapse 
process. As a wave travels outward, its (approximately 
constant) physical width corresponds to fewer coordinate 
grid spacings. Eventually, the wave fails to be resolved, 
and since the system is mixed hyperbolic-elliptic, this 
failure of resolution in one part of the grid can affect the 
entire grid. This difficulty means that for a given spatial 
resolution, there is only a certain amount of time that 
the evolution can be run and still give reliable results. 
This places a limit on the type of problems that can be 
treated using this method. 

The variables (S, K z r , U, W) are evolved using a 3 step 
iterative Crank-Nicholson (ICN) algorithm. At each step 
of the ICN process, the elliptic equations for (ip, a, (3 r , (3 Z ) 
are solved using the conjugate gradient method with Neu- 
mann preconditioning. 

Though we solve for ip using the Hamiltonian con- 
straint, we note that ip also satisfies an evolution equa- 
tion. From equation (0) we have 



d t ip = p z d z ip + (3 r d r ip + ip 



r a 

l + - Q (U + 2rW) 



(15) 



We use equation ( |l5| ) as a check on the reliability of the 
code. Specifically, we compute a second value of ip by 
starting with the initial data and then evolving using 
equation (|l||). We then check that the two values of ip 
agree well at r = z = 0. We run the code only as long as 
this good agreement persists. 

The convergence of the code is tested using the mo- 
memtum constraint p a = D h K a \ ) (Since K = 0). The 
Hamiltonian constraint is solved for ip and therefore can- 
not be used as an independent test. In general, the meth- 
ods we use would lead one to expect second order conver- 
gence. However, the boundary conditions at infinity for 
the elliptic equations may be only first order accurate. 
Figure 1 shows the L 2 norm of p z plotted as a function 
of time. For these runs, a = 4 and ay = a z = 1. The 
dashed line corresponds to 42 gridpoints in the r direc- 
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tion and 42 gridpoints in the z direction, while the solid 
line corresponds to a run with twice the resolution. The 
results indicate second order convergence. 
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FIG. 1. IpJ plotted vs. t for two different resolutions 



III. RESULTS 

The code was run in double precision on Dec alpha 
workstations and on the NCSA Origin 2000. Unless oth- 
erwise stated, all simulations were run with 42 gridpoints 
in the r direction and 42 gridpoints in the z direction. 

The initial data that we use, and therefore the space- 
time that it evolves to, has reflection symmetry about 
the z = plane. Though our algorithm does not require 
this extra symmetry, since the spacetimes we treat have 
this reflection symmetry, we save computational time by 
running the simulations in the range < z < ir/2. The 
reflection symmetry and axisymmetry give spacetime a 
preferred world line: the origin (r = z = 0). In gravi- 
tational collapse in maximal slicing, it is expected that 
the lapse approaches zero in the region of strong gravity. 
Therefore, a useful quantity to plot is In cxq as a function 
of time, where cto is the lapse at the origin. Note that this 
quantity can be directly compared between two different 
codes, provided that both use maximal slicing. 

Our initial data has three parameters: the amplitude a 
and the two widths a r and a z . However, it turns out that 
the effective parameter space is two dimensional. Under 
the transformation ay — > cay, a 2 — > ca z , a — > a/c 2 for a 
positive constant c, the spacetime is changed only by an 
overall constant scale. We can think of the two remaining 
parameters as being the strength and shape of the wave, 
and we want to find the dependence of the collapse on 
these two parameters. Figure 2 shows the dependence of 
the collapse on the strength of the wave. Here, three sim- 
ulations are run, each with ay = a z = 1. The parameter 
a is 4 for the solid line, 5 for the dashed line and 6 for the 
dot-dashed line. The a — A and a = 6 spacetimes have 



been studied using a 3+1 code in reference |L7|]. Our re- 
sults are in good agreement with theirs. For a = 4, the 
lapse, after an initial collapse to small values, appears 
to be returning to I. That is, we expect the Brill wave 
to disperse. For a — 6, the lapse continues to collapse, 
and we expect a black hole to form. Somewhere between 
a = 4 and a = 6 is an amplitude that leads to critical col- 
lapse; however our method does not have the resolution 
needed to study this process. The reason for this is that 
to study critical collapse, one must evolve long enough 
to distinguish a spacetime slightly below the black hole 
threshold from one that is slightly above. In this amount 
of time, the waves traveling towards the outer boundary 
will, in general fail to be resolved, leading to a general 
lack of accuracy in the results of the evolution. 
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FIG. 2. lnao plotted vs. t for three different amplitudes: 
a — 4 (solid line), a — 5 (dashed line) and a = 6 (dot-dashed 
line) 



The issue of the dependence of the collapse process 
on the shape is somewhat subtle, due to the question of 
what is to be held constant as the shape varies. The sim- 
plest thing to do is to hold the amplitude a and the sum 
ay + u z constant while varying a r — a z . Results of these 
collapse simulations are shown in figure 3. Here, a = 4 
and a r + <r z = 2 for all simulations. We will call initial 
data "spherical" if <r r — a 2 , "prolate" if a r < a z and 
"oblate" if ay > a z . (This terminology is somewhat mis- 
leading, since it is not the wave itself but the initial S/ r 
that is spherical, prolate or oblate). Figure 3 contains a 
spherical collapse (solid line), a prolate collapse (dashed 
line) with a r — a z = —0.4 and an oblate collapse (dot- 
dashed line) with a r — <r z = 0.4. Here the shape seems to 
have a large influence on the collapse, with even a small 
amount of oblateness producing collapse, while a small 
amount of prolateness hastens dispersion. 

However, it is not clear how to interpret this result, 
since for a fixed a and a r + a z , it may be that there is 
a dependence of the "strength" of the gravitational wave 
on o r — g z . To see this, we examine how the ADM mass 



4 




1 2 3 4 5 

t 

FIG. 3. lnao plotted vs. t for a = 4 and three different 
shapes: spherical (solid line) , prolate (dashed line) and oblate 
(dot-dashed line) 



depends on the degree of prolateness or oblateness at 
fixed a and a r + cr z . These results are shown in figure 4. 
Here, we have ay + a z = 2 and a = 4. The ADM mass M 
is plotted as a function of ay /ay. Since, at a fixed a, M 
is an increasing function of a r /a Zl this is (at least part 
of) the reason why at fixed a oblateness causes collapse. 

2 I — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — 
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FIG. 4. M plotted vs. a r /cr z for a = 4 and oy + a z = 2 

One might therefore instead look at the influence of 
shape at fixed oy + a z and fixed M. The results of this 
comparison are shown in figure 5. Here oy + a z = 2 and 
M = 0.46 (the value of M corresponding to a = 4 and 
ay = ay = 1). The solid line is a spherical collapse, while 
the dashed line is a prolate collapse with ay — ay = —0.4 
(corresponding to a = 6.1) and the dot-dashed line is an 
oblate collapse with ay — ay = 0.4 (corresponding to a = 



2.7). Here the shape has some influence on the details of 
the collapse process; but this much change in the shape 
does not seem to have a particular tendency either to 
promote or to inhibit collapse. Figure 6 shows the same 
sort of plot, but with more distortion in the shape. Again 
M = 0.46 and oy + oy = 2. However, here the prolate 
collapse (dashed line) has a r — o 2 = — 1 [a — 14) and 
the oblate collapse (dot-dashed line) has o r — o 2 = 1 
(a = 1.4) in addition to the spherical collapse (solid line). 
Here it seems that at constant ADM mass there is a slight 
tendency of oblateness to induce dispersion. 
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FIG. 5. lnao plotted vs. t for M = 0.46 and three different 
shapes: spherical (solid line) , prolate (dashed line) and oblate 
(dot-dashed line) 
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FIG. 6. lnao plotted vs. t for M = 0.46 and three differ- 
ent shapes: spherical (solid line), prolate (dashed line) and 
oblate (dot-dashed line). Here the amount of oblateness or 
prolateness is greater than that in figure 4 

We now consider more strongly gravitating initial data 
and follow its evolution until the formation of an appar- 
ent horizon, while considering some of the properties of 
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the collapse. This simulation has M = 2, ay = 1 and 
a z = 1 (a = 8.5) and is run with 82 gridpoints in the 
r direction and 82 gridpoints in the z direction. Figure 
7 shows lnao as a function of time. The lapse collapses 
throughout the evolution. At each time step we calculate 
the Riemann invariant I — R abcd R a bcd/16 and find the 
maximum of its absolute value Z max as well as the spatial 
position where /| = / max . Figure 8 shows ln/ max plotted 
as a function of time. Here we see that after an initial 
increase, In J max decreases during the rest of the evolu- 
tion. At all times during the evolution the place where 
\I\ = W is the origin. 

To monitor the approach to apparent horizon forma- 
tion, we consider how the horizon is found. On a max- 
imal slice in an axisymmetric spacetime, the horizon is 
given by a curve z = f(r) where the function / satisfies 
a differential equation that is integrated from the axis 
to z = 0. Let 9 denote the angle at which the curve 
meets the z = plane. Then the curve is a horizon only 
if 9 = 7r/2. The horizon finding subroutine integrates 
the differential equation for / starting at each point on 
the axis, finds the angle 9 for each curve and then finds 
^max, the maximum value over all curves of this angle. If 
$max < 7r/2 then there is no horizon at this time. Figure 
9 shows 28 max /-jv plotted as a function of time. Note that 
Omax increases throughout the collapse process. Note also 
that even before the actual horizon formation, one can 
tell that a horizon is about to form by noticing that # m ax 
is approaching n/2. The horizon forms at t — 3.9 with 
area A = 165 = 0.82(16ttM 2 ). 
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FIG. 7. lnao plotted vs. t for M = 2 and er r = a z = 1 



We now turn to a highly prolate collapse: the evolu- 
tion of one of the initial data sets of reference . Here, 
M = 2, cr 2 = 1.6 and a r = 0.128 (a = 325). This simula- 
tion was run with 162 grid points in the r direction and 
42 grid points in the z direction. Ideally, we would like to 
follow the collapse until either a horizon forms or a cur- 
vature scalar blows up. Unfortunately, we are not able to 
follow the evolution that long; so we evolve for a some- 
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FIG. 9. 29vwc/ir plotted vs. t for M = 2 and ay = a z = 1 



what shorter time and attempt to discern trends in the 
evolution. Figure 10 shows In ao plotted as a function of 
time. Here we see the usual collapse of the lapse through- 
out the evolution. Figure 11 shows ln/ max plotted as a 
function of time. Here we see that the maximum of the 
Riemann invariant decreases as the collapse proceeds. In 
the initial data, the spatial location where |/| = / ma x is 
on the axis at z = 1.04. As the evolution proceeds, this 
spatial location remains on the axis, but moves towards 
the origin, reaching the origin at t = 0.55 and then re- 
maining at the origin for the rest of the evolution. To 
discern a trend in the approach to apparent horizon for- 
mation, we plot (Figure 12) 2(9 max /7r as a function of time 
for this simulation. Note that this quantity is increasing. 

The trends of this evolution are that (i) 7 max decreases, 
(ii) the spatial position where |/| = / max moves to the 
origin and (iii) # max increases. If these trends continue, 
then this spacetime will form a black hole rather than 
a naked singularity. Thus (in this example at least) it 
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FIG. 10. ln« plotted vs. t for M = 2, oy = 0.128 and 
o z = 1.6 



FIG. 12. 20 max /7r plotted vs. t for M = 2, oy = 0.128 and 
o z = 1.6 
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FIG. 11. ln/ max plotted vs. t for M = 2, oy = 0.128 and 
a 2 = 1.6 



seems that Brill waves behave differently from collision- 
less matter with even highly prolate initial configurations 
forming black holes. This conclusion is not firm for two 
reasons: (i) we have only followed the evolution for a 
certain amount of time, and the trends that we have ob- 
served in this part of the evolution could reverse in later 
parts, (ii) we have only evolved a certain, highly prolate 
initial configuration. It is possible that much more pro- 
late initial configurations behave differently. Both these 
issues need further study. 

IV. CONCLUSION 

We have presented (i) a new method of curing axis in- 
stabilities, (ii) a study of the dependence of Brill wave 
collapse on the shape of the wave including (iii) a pre- 
liminary examination of the issue of cosmic censorship in 
the collapse of highly prolate Brill waves. A convincing 



demonstration that the axis really is stable would require 
evolution for long times. However, our outer boundary 
condition entails a loss of resolution that increases with 
time and this prevents us from running the code for a 
long time, even for weak waves. (Note: one might also 
expect that our sort of outer boundary condition would 
result in spurious reflected waves. While we cannot run 
our code long enough to check this, we have used our 
outer boundary condition to perform a numerical simu- 
lation of the linear wave equation; and there we find very 
little reflection). 

A way to test the regular axis method independently of 
the question of outer boundary conditions is to treat the 
case of closed cosmologies. One of us has used the regular 
axis method to treat Gowdy spacetimes on S 2 x S 1 x R. 
fLU Here the evolution proceeds for long times and there 
is no axis instability. Though the Gowdy spacetimes have 
two Killing fields, one can easily apply the regular axis 
method to closed axisymmetric cosmologies. 

To do a more thorough study of Brill wave collapse, we 
require a better way of treating the outer boundary. One 
way to do this is to put the boundary at a finite distance 
but replace the ADM system by one in which a simple 
outgoing wave boundary condition does not cause nu- 
merical instability. The method of Baumgarte, Shapiro, 
Shibata and Nakamura jn|^0| seems to have the appro- 
priate properties. Another method would use harmonic 
coordinates, which makes Einstein's equation similar to 
the wave equation and should be stable with an outgoing 
wave boundary condition. Alternatively, one could keep 
spatial infinity (or null infinity) on the computational 
grid, but replace the ADM equations with the conformal 
equations of Friedrich. |2l],|22|,|23| We are presently exam- 
ining these alternatives to see which is likely to work best 
for Brill wave simulations. 

Our preliminary results on Brill wave collapse indicate 
that even highly prolate Brill waves do not evolve to form 
naked singularities. In the evolution of highly prolate ini- 
tial data, gravitational collapse will tend to increase the 
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Riemann curvature, while the tendency of the wave to 
spread out will lessen curvature. We have shown that in 
the initial stage of the evolution of a highly prolate wave, 
the second effect is the more important one. That is, the 
highly prolate wave tends to become less prolate rather 
than more as it evolves. In this way, Brill waves behave 
differently from the collisionless matter studied in refer- 
ence [[To). To put our tentative conclusions on a firmer 
footing, we need to study more extreme configurations 
and evolve them for a longer time. 
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